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ABSTRACT 

We study the effects of a large-scale, ordered magnetic field in protoplanetary disks on Type I 
planet migration using a linear perturbation analysis in the ideal-MHD limit. We focus on wind¬ 
driving disks, in which a magnetic torque oc Bo^dBoip/dz (where i?oz and Bq^ are the equilibrium 
vertical and azimuthal field components) induces vertical angular momentum transport. We derive 
the governing differential equation for the disk response and identify its resonances and turning points. 
For a disk containing a slightly subthermal, pure-Boz field, the total 3D torque is close to its value 
in the 2D limit but rema ins lower than th e hydrodynamic torque. In contrast with the 2D pure-i?ov 3 
field model considered by Terquem (20031, inward migration is not reduced in this case when the field 
amplitude decreases with radius. The presence of a subdominant Bg^p component whose amplitude 
increases from zero at z = 0 has little effect on the torque when acting alone, but in conjunction with 
a Bqz component it gives rise to a strong torque that speeds up the inward migration by a factor 
> 200. This factor could, however, be reduced in a real disk by dissipation and magnetic diffusivity 
effects. Unlike all previously studied disk migration models, in the Bqz + dBQ^jdz case the dominant 
contributions to the torque add with the same sign from the two sides of the planet. We attribute 
this behavior to a new mode of interaction wherein a planet moves inward by plugging into the disk’s 
underlying angular momentum transport mechanism. 

Subject headings: accretion, accretion disks - MHD - planet-disk interactions - protoplanetary disks 


1. INTRODUCTION 

Type I migration is the radial drift induced in a planet 
that is embedded in a protoplanetary disk by its linear 
gravitational interaction with the surrounding gas. This 
process, which is applicable to low-mass planets (typi¬ 
cally up to Neptune’s mass), is of great relevance to the¬ 
ories of plan e t formation and evolu t ion (see, e.g.,|Lubow 
&: Ida 2010[ [Kley fc Nelson||201^ Baruteau fc Mas^ 
2013 aM Baruteau|2013|for recent reviews). Early treat- 
ments of this interaction focused on isothermal disks with 
a smooth density distribution and inferred rapid inward 
migration. The subsequent discovery of “hot Jupiters,” 
giant planets orbiting within ~ 10 stellar radii of their 
host stars, has exacerbated the conundrum elicited by 
this finding. These large, mostly gaseous, planets must 
have formed much farther out in their natal protoplan¬ 
etary disks, where, according to the currently favored 
core-accretion model, they needed to assemble a large- 
enough (> a few times M®) rocky core before a fast gas 
accretion phase (which created their envelopes) was trig¬ 
gered. However, the rapid inward drift predicted by the 
early models of Type I migration would prevent the pos¬ 
tulated embryonic cores from lingering long enough at 
their formation sites for consistency with the observed 


distribution of hot Jupiters (e.g., Ida & Lin 

2008 

I- 

In the companion paper (Uribe et al. 

W 


hereafter 


Paper I) we briefly review some of the more recent at¬ 
tempts to study Type I migration in the context of more 
realistic protoplanet ary disk models, and we note a few 
of the promising results that have already been obtained. 
These include the effects of a disk magnetic field, which 

^ present address: Astronomy Department, Adler Planetar¬ 
ium, Chicago, IL 60605, USA 


could be either a small-scale, disordered field associated 
with magnetohydrodynamic (MHD) turbulence (in par¬ 
ticular, the turbulence induced by the magnetorotational 
instability) or a large-scale, ordered field that permeates 
the disk. These two field configurations likely play a cen¬ 
tral role in, respectively, the radial and vertical an gular 


bus 

2011 


Konigl & Salmeron 

2011 

). The effects of a 

sma 

ll-sca 

e 

, disordered tield on planet migration include 


the stochastic motions induced in the smallest planets 
by MHD turbulence, the effective viscosity that prevents 
saturation of the corotation torque arising from the co¬ 
orbital region (which can be the dominant component of 
the torque on a planet undergoing Type I migration in 
a nonisothermal disk), and the magnitude of the mass 
threshold for opening a gap in the disk (i.e., for the on¬ 
set of Type H migration), which depends on the value 
of the effective viscosity. A novel feature of the pres¬ 
ence of a large-scale, ordered field is the appearance of 
new types of waves that propagate in the disk (specif¬ 
ically, slow- and fast-magnetosonic [SMS and FMS, re¬ 
spectively] , as well as Alfven), which supplant the sound 
waves that characterize a hydrodynamic (HD) disk. The 
MHD waves modify the way in which the disk responds 
to the gravitational perturbation induced by the planet, 
and, in turn, the way in which the disk acts back on 
the planet to cause it to migrate. The studies conducted 
so far have demonstrated that a purely azimuthal field 
in 2D can slow down, and even reverse, inward migra¬ 
tion if the magnetic field pressure is not much smaller 
than the thermal pressure and i f the field amplitude de¬ 
creases with radius fast enough ( Terquem|2003 Fromang 


et al. 2005). It was also shown that, while a strong az- 
imuthal held suppresses the corotation torque, a weaker 
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one modifies it in a w ay that can potenti ally lead to out- 


ward migration (e.g., 

Gullet et al. 

2013). In the case of 

a purely vertical field it was found by 

Muto et al. I 

20081 


that the only effect in 2D is the replacement ot the HD 
sound waves by the MHD FMS waves, with a resulting 
reduction in the magnitude of the torque, but that, in 
3D, SMS and Alfven waves also play a role (although 
these authors could not determine the net torque since 
their calculations were carried out in the shearing-sheet 
approximation). 

As we observed in Paper I, a field configuration that 
is either purely azimuthal (and uniform with height) or 
purely vertical, as assumed in the aforementioned papers, 
is not representative of real protoplanetary disks. Such 
disks are likely permeated by open field lines that cor¬ 
respond to the interstellar field that originally threaded 
the parent molecular cloud cores and was dragged in¬ 
ward when the cores underwent gravitational collapse 
and formed the central stars and their associated circum- 
stellar disks. The radially advected field lines assume a 
pinched (or hourglass) morphology and, in addition, are 
twisted by the differential rotation in the disk. In gen¬ 
eral, the field would thus have vertical, radial, and az¬ 
imuthal components, which, near the disk surfaces, could 
have comparable magnitudes. This morphology is con¬ 
ducive to the formation of centrifugally driven outflows 
that can efficiently transport angular momentum away 
from the disk. This picture is supported by observations 


Konigl & Salmeron 

2011 

as Bans & Konigl 

2012 


and references 


toplanetary disk, the magnetic field morphology in the 
giant-planet formation region can be expected to main¬ 
tain an equilibrium structure on a timescale that greatly 
exceeds the local rotation period because the bulk of the 
gas at that distance (a few AU) is generally weakly ion¬ 
ized and therefore magnetically diffusive. In the simplest 
representation of such an equilibrium open-field config¬ 
uration, the field has an even symmetry about the mid¬ 
plane (z = 0 in cylindrical coordinates {r^ip,z})^ with 
the equilibrium (subscript ‘0’) radial (i?or) and azimuthal 
(Boip) field components changing sign there (but with the 
vertical component Bq^ remaining nonzero). 

The aim of Paper I and this paper is to conduct a sys¬ 
tematic investigation of the effects of a multi-component, 
ordered, large-scale magnetic field on Type I planet mi¬ 
gration. The ultimate goal of this study is to model a 
planet that is embedded in a wind-driving disk, which, as 
the discussion above indicates, involves all three spatial 
components of the field and must, for self-consistency, be 
treated within the framework of nonideal MHD. Given 
the complexity of this problem, its full treatment is de¬ 
ferred to a future publication, and the present study con¬ 
centrates on simpler field configurations that can be in¬ 
vestigated using ideal MHD. We employ the complemen¬ 
tary approaches of a linear perturbation analysis (the 
focus of this paper) and of numerical simulations (the 
main focus of Paper I). The linearization procedure is 
formulated in full generality in Section of this paper, 
the numerical procedure employed in calculating the net 
torque on the planet is described in Section]^ and semi- 
analytic solutions in 2 and 3 dimensions for three repre¬ 
sentative field configurations are presented in Section]^ 


The most general field configuration that can be treated 
self-consistentlv within the framework of ideal MHD in¬ 
volves a field that has both vertical and azimuthal com¬ 
ponents that can varv with r but not with z. This 
case is discussed extensively in Paper I (which, in Sec¬ 
tion 1.2, also presents analytic results obtained from our 
perturbation analysis for this case)jj The 3 examples 
considered in this paper are a pure-B^ field, a pure-H,^ 
field whose amplitude increases with height (from zero 
at z = 0), and the two combined. The first two con¬ 
figurations are treated numerically in Paper I. The third 
configuration, which induces vertical angular momentum 
transport, mimics the situation in wind-driving disks. It 
is, however, not entirely self-consistent in the context of 
ideal MHD since the radial velocity that results from the 
loss of angular momentum would give rise to a radial 
magnetic field component that, in the absence of mag¬ 
netic dilfusivity, would be rapidly sheared out of equi¬ 
librium by the differential rotation in the disk. This ex¬ 
ample is, nevertheless, included here in the hope that it 
could provide useful insights into the more general prob- 
lemld^Finally, Sections^(complemented by an appendix) 
andprovide a discussion and a recapitulation, respec¬ 
tively, of the results presented in this paper. 

2. LINEAR ANALYSIS 

2.1. Basic Equations 

The dynamics of a magnetized protoplanetary disk is 
governed by the momentum conservation, mass conser¬ 
vation, and induction equations, given respectively by 


U\ 

p — -|-(v-V)v =-VP-1-F - pVt/f , 

(1) 

| + v-(rv) = o, 

(2) 

and 


-^-hVx(vxB)=0, 

(3) 

where p is the mass density, ijj = GMp/\r — rp the grav¬ 
itational potential (with Mp and rp being the planet’s 
mass and orbital radius, respectively and G the gravi¬ 
tational constant), v the velocity field, B the magnetic 
field, and F the Lorentz force density 

F = (V X B) X B/p 

(4) 


(with = 47r). The field is required to satisfy the 
solenoidal condition V • B = 0. 

The induction equation was written in its ideal-MHD 
form, which deserves a clarification. As already noted in 
Section IB the giant-planet formation zone in protoplan¬ 
etary disks is expected to be weakly ionized and hence 
magnetically diffusive. This observation applies in par¬ 
ticular to the disk midplane, where the planets reside, 
since this region is most strongly shielded from the main 
ionization sources of the disk (external cosmic rays, X- 
rays, and UV radiation). As a result, the ideal-MHD 

^ We henceforth preface section, equation, and figure numbers 
in Paper I by the numeral ‘L. 

® A similar motivati on u nderlies the discussion of a vertical gra¬ 
dient of Br in Section |2.6[ 
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limit can be employed only under certain constraints. 
Since the focus of the current study is the effect of the 
magnetic held on the torque that the disk exerts on the 
planet (which arises from the nonaxisymmetric pertur¬ 
bations that the planet induces in the disk), it seems 
reasonable to require that the effective coupling time 
between the neutral disk component (which dominates 
the gravitational interaction with the planet) and the lo¬ 
cal magnetic held be shorter than the Keplerian rotation 
period (since any nonaxisymmetric perturbations would 
be washed out over timescales much larger than the in¬ 
verse of the Keplerian angular velocity Dk). The ratio of 
these two timescales is quantihed by the Elsasser number 
A = where v\ = /^ip is the square of the 

Alfven speed and r]± is t he magnetic diffusiv ity perpen¬ 


dicular to the held (e.g., Konigl et al. 


2010). It can be 


expected that a magnetized disk would influence planet 
migration in an ideal-MHD-like way if A 3> 1 at the 
midplane. Perhaps not surprisingly, this is also the min¬ 
imum coupling condition for sustaining angular momen¬ 
tum transport through either MRI -induced turbulence 


or a centrifugally driven wind (e.g., Konigl & Salmeron 


20111. It is conceivable that the diffusivity ot a real disk 


could give rise to interesting effects that cannot be cap¬ 
tured with the current formulation, but it nevertheless 
seems worthwhile to study the simpler ideal-MHD limit 
hrst. 

Throughout this work we use an inertial, nonrotat¬ 
ing, cylindrical coordinate system (r, tp, z) centered on 
the star. The adopted held geometry and equilibrium 
disk conditions are described below. We also adopt an 
isothermal equation of state, P = c^p, where c is the 
isothermal speed of sound. Although our focus is on 
the disk midplane, we allow for vertical wave propaga¬ 
tion, which enta ils taking account of the 3D structur e 
of the disk (e.g., Tanaka et al.| 2002 Muto et al. 2008). 
We simplify the treatment by averaging over the ettec- 
tive thickness 2h of the disk (assumed to be <C r at the 
location of the planet) while imposing an even held sym¬ 
metry. In this approach, the density p is replaced by the 
surface density E, although we will also be referring to 
the verticall y averaged equilibr ium density Pav = Eo/2/i. 

Following Muto et al. (2008), we further simplify our 
treatment by neglecting the effects of the vertical com¬ 
ponent of the stellar gravitational held on the disk struc¬ 
ture. 

In the remainder of this paper we normalize the spatial 
variables r, z, and h by the orbital radius of the planet 
Cp, frequencies by the planet’s orbital frequency Dp, time 
by Dp and velocities by rpDp. 


2.2. Equilibrium Conditions 

We adopt a simplihed equilibrium magnetic held con- 
hguration intended to capture the basic geometry of an 
even-symmetry open magnetic held that threads a geo¬ 
metrically thin disk: 


Bo = {Bor{z),Bo^{z),Bo^{r)} . (5) 

The geometrical thinness h <C I implies, by the 
solenoidal condition, that Bqz inside the disk is very 
nearly constant with height. Under the assumed even 
symmetry, the radial and azimuthal held components 
vanish at z = 0 and increase in magnitude on going away 


from the midplane (with opposite signs in the two hemi¬ 
spheres). Based on a simplihed treatmen t of a wind¬ 
driving disk (e.g., Wardle & Konigl 1993), we approx¬ 
imate Boip(z) = i3^o^7”wTiere~5(^is a constant, and 
similarly for Bor{z). Since the dominant contribution to 
B,^ in such disks is the shearing of Br by the differen¬ 
tially rotating gas, B^o generally has the opposite sign of 
Bro- We henceforth concentrate on the case where Bqz 
and Bor are > 0 and i?o<^ is < On We adopt a power- 
law form Boz(r) = (with B^o a constant) for the 

radial dependence of the vertical held component, which 
is motivated by self-similar models of wind-driving disks 


(e.g., Blandford & Payne 1982 Contopoulos & Lovelace 


19941. 'I'he other equilibrium magnetic held components 

would in general also vary with r, but we neglect this 
dependence in this paper. (It is, however, incorporated 
into the functional form we adopt for Bq^ in Paper I.) 

Using the above held morphology in Equation ([4]), we 
obtain for the vertically averaged equilibrium Lorentz 
force 


p{F) = (^BroB^or^^ - qzBiy‘1^-^ - B\ 


2 


B^oBzo'f'^’' + B^oBrQ— , 0 


( 6 ) 


where the triangular brackets () denote a vertically aver¬ 
aged (between —h and h) quantity. The hrst two terms 
in the radial component of this expression represent the 
magnetic tension and pressure forces, respectively, which 
act to reduce the inward pull of gravity. (The third ra¬ 
dial term is generally much smaller due to the h? factor.) 
As we noted in Section the strong shearing of a radial 
held component cannot oe self-consistently incorporated 
into an equilibrium ideal-MHD model, so, in the applica¬ 
tions given in this paper and in Paper I, we set i?or = 0. 
This eliminates the radial tension, which would typically 
dominate the radial Lorentz force in a real wind-driving 
disk. The hrst (and dominant) term in the azimuthal 
component of Equation)]^ represents the braking torque 
acting on the disk. When substituted into the (p compo¬ 
nent of Equation ([^ (which describes the conservation 
of angular momentum), we obtain an expression for the 
induced accretion velocity vor. 


_ 2D(F^) 

^Or — 9 5 

Pavl^^ 


(7) 


where D is the angular velocity and k is the epicyclic 
frequency (k^ = 

Equation ([^ highlights the fact that a disk threaded 
by a vertical magnetic held and possessing a nonzero ver¬ 
tical gradient of the azimuthal magnetic held would ex¬ 
perience vertical angular momentum transport. This is 
a key property of disks that drive centrifugal outhows. 
Because of the strong vertical gradients that character¬ 
ize geometrically thin disks, this transport could be very 


^ Note that our ansatz implies Bipo = dBipjdz, so the constant 
B^O can be used interchangeably with dBtp/dz when referring to 
the effect of a vertical gradient of the azimuthal field component. 
We distinguish this effect from that of an azimuthal field that is 
finite at the midplane and independent of z, which is considered in 
Paper I. 
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efficient. In fact, when such a wind is launched, the 
resulting accretion speed (assuming adequate magnetic 
coupling, A >> 1, in the disk) is of the order of c (e.g., 
Salmeron et al. 2007), which is much higher than the 
speeds that arise trorn radial angular momentum trans¬ 
port by a small-scale, disordered field. As was noted in 
Section the field-line bending induced by the inflow¬ 
ing gas would give rise to a radial field component that, 
in turn, would lead to a rapidly growing azimuthal field 
component and a departure from equilibrium unless mag¬ 
netic diffusivity were taken into account. Since we intend 
to include nonideal-MHD effects in a future work, we re¬ 
tain the VQr term in the basic formulation presented in 
this paper and assume an equilibrium velocity field of the 
form V = (uor,rU,0) in the ensuing perturbation analy¬ 
sis. However, the semianalytic results that we obtain are 
derived in the limit Vq^ = 0, which should be an adequate 
approximation given that c/rfl (which is of the order of 
h) is < 

2.3. The Disk Response 

To explore how our equilibrium disk model responds 
to the perturbing potential of the p lanet, ip' (norma lized 
by r-pUp), we follow the method of Terquem (2003) and 
earlier analyses by considering small perturbations and 
linearizing the basic equations. We allow the planet to 
perturb the disk in all 3 directions, and Fourier trans¬ 
form in time as well as in space along the tp and z di¬ 
rections by considering all Fourier quantities to be of the 
form A(^ with the frequency taken 

to be a; = mflp. As in Terquem (2003), we denote 
all perturbed quantities with a prirnel The planet’s 
potential itself can be expanded in a Fourier series, 
•0' = X]m=o '^'m cos mcf), where (j) = ip — flpt, and ex¬ 
pressed in te rms of the generalized Laplace coefficients 
bjUa) (e.g., 


^r/2(«) 

Terquem]|2003‘f T 


Ward 


1989 


Korycansky & Pollack 


1993 


, _ I r > 1 


( 8 ) 


where M* is the stellar mass (assumed to be ^ Mp), the 
coefficient 6^ is equal to 0.5 for m = 0 and to 1.0 for all 
the other (positive integer) values of m. 


'' 1/2 


2 r 
(r) = - 

^ Jo 


cos m(j) d(p 


lo (<7^— 2acos((>)i/2 ’ 

and where, for r > 1, a = 1/r, q = 1, and = 1 -I- e^; 
whereas for r < 1, a = r, p = 1, and = 1 -|- (with 
e, a real number of magnitude <C 1, being the potential 
softening parameter; see Equation (1.20)). 

Linearizing equations ([^ and ([^ yields 


imav'j. — 2v'^fl + 


dK , / ^'COr 

Vor^^ + Vr 


, dvo 


dr 

2h 


dr 

W' 


dz 


-Yr^^' + W') + ^y^{F,) + {F',)). (10) 

^ A possible exception to this con clusi on is our treatment of 
the case 0 in Section |4.4| which, as was already 

pointed out in Section^ does not handle tne equilibrium magnetic 
field configuration in a self-consistent manner. 


imav, 


fdv' 




+ F {F')) , (11) 

Zjq C 


dv' 


ip! 

{F') 


imaov' + = —ikzW — ik^ip' + 2,h 

dr So 


( 12 ) 


and 


imaW' —Id, ,, 
- 2 — = 

rLo dr ^ ‘ 


1 d 


- —< - ~^p'Wr ’ d3) 


where W = E'c^/Eq is the enthalpy perturbation, a = 
U — Up (= U — 1 in normalized units), and where we 
suppressed the Fourier labels (m and kz) on the primed 
quantities to simplify the presentation. To calcula te t he 
pert urbed Lorentz force components in Equations (10) - 
(12), we use 


pF' = (V X B') X Bo -f (V X Bo) x B' (14) 

and the foll owing expression from [Chandrasekhar ^ 
Fermi (1953), 


B' = V X (I X Bo) , 


(15) 


which relates the perturbed magnetic field to the La- 
grangian displacement This relation expresses the 
ideal-MHD condition of the magnetic field lines being 
“frozen” into the matter, and can be derived from the 
induction equation in t he l imit of a negligible equilib¬ 
rium velocity. Equation (15) does not strictly hold when 
Vor is finite, but its use in the present analysis is justified 
by the fact that the results we derive are obtained in the 
limit VQr = 0. 

Note that the vertical averaging in the perturbation 
equations is done after all the spatial derivatives are eval¬ 
uated. Thus, for example, the perturbed vertical mo¬ 
mentum equation contains a term oc dB^^p/dz (if it is 
nonzero) even though the equilibrium vertical magnetic 
force vanishes in our model (see Equation ([^). It is also 
worth noting that, in general, the perturbation equations 
do not contain terms involving the spatial derivative of 
the equilibrium density if the temperature is assumed 
to be a s patial constant (as is done in th is paper, follow - 
ing, e.g., Korycansky & Pollack 1993 and Terquem|2003 ). 
Consequently, the linearization results are not sensitive 
to our neglect of a vertical density gradient (induced by 
either tidal gravity or a vertical magnetic pressure gra¬ 
dient) in the equilibrium momentum equation. 

The Lagrangian displacement is related to the Eulerian 
velocity perturbations through 


'' = a+<''» 


V)4-(e-V)vo 


(16) 


(e.g., Zhang & Lovelace 


now on to suppress the 
quantities. 


2005|. Thus, continuing from 
i'ourier labels on the perturbed 


/ ■ p dVor p 

= ima^r + '^Or-^ -’ 

V = ima^p + VQr -p. - 


dr 


dr 


(17) 

(18) 
















































Type I Planet migration in a Magnetized Disk. 11. 


= ima^z + VQr -T- 
or 


(19) 


Substituting these expressions into Equations (101 - 
and carrying out the necessary algebra, one obtains a 
second-order differential equation for of the form 

= Si{r)^ + Soir)i^' , 
(20) 

which can be integrated numerically (see Section 
The derived solution for can then be plugged bacK 
into the linearized equations of motion and continuity 
to obtain W and to evaluate the torque exerted by the 
planet on the disk. 

2 . 4 . Resonances and Turning Points 


The solution of Equation (20) becomes singular at the 
locations where the leading-order coefficient, A 2 (r), goes 
to zero. These locations correspond to the resonances of 
the differential equation and are manifested by a diver¬ 
gence of the density perturbation. Their presence also 
signals the appearance of distinct regions of wave prop¬ 
agation in the disk. The torque exerted in the regions 
surrounding the resonances can potentially dominate the 
response of the disk and thus determine the rate and di¬ 
rection of planet migration. It is therefore important to 
calculate the resonant locations and examine their con¬ 
tribution to the integrated torque. 


In general, the differential equation (20) describes wave 


propagation in the disk, and we can rewrite a homoge¬ 
neous version of it in terms of just th e second- and zer oth- 
order terms. We do this, following Terquem (20031, by 


defining a new dependent 


variable, which yields 

d^Xr 


dj.2 


+ KXr = 0 , 


where 


A2 4[A2J 2dr[A2 


( 21 ) 


( 22 ) 


When K is real. Equation (21) has wave-like solutions if 


K > 0 and evanescent solutions if if < 0. More gener¬ 
ally, when K is complex (with real and imaginary parts 
given by and 3{itr}, respectively), this equation 

always has wave-like solutions, but whether or not the 
waves propagate is determined by the r nag nitude and 
signs of and (see Appendix [^for details). 

In our model, K is complex for disks tnat contain a 
vertically variable equilibrium field component (nonzero 
dBip/dz or dBr/dz), and real for vertically constant field 
configurations. The locations where the solution changes 
its character from wave propagation to evanescence (i.e., 
where K = 0) are known as the turning points of the 
differential equation. Since it is the density wave propa¬ 
gation from the turning points that is responsible for the 
exchange of angular momentum with the planet (e.g., 

® The functions An. A i. A 2 , Sq, and Si that appear as co¬ 
efficients in Equation ( |20[ | are provided as online supplementary 
material in the form of both a Mathematica notebook and a PDF 
document. 


Goldreich fc Tremaine|1979 |, finding the locations of the 
turning points makes it possible to identify the regions 
of the disk that exert a back torque on the planet. In 
the next subsection we evaluate the locations of the reso- 
nanc es a nd turning points of a simplified version of Equa¬ 
tion (201 and show how their properties depend on the 


magnetic field configuration. 


2 . 5 . Two- and Three-Dimensional Modes 


With a nonzero vor, Equations (10) - (13) all contain 
both zeroth-order and derivative terms of the perturbed 
variables, which means that it is impossible to eliminate 
all variabl es be sides from the equations. As discussed 
in Section |2.2[ vgr is an effectively small velocity, so in 
this initial study we drop all the t erm s in which it ap¬ 
pears in the equations. Equation (15) then manifestly 
satisfies the linearized induction equation, and the (ver¬ 
tically averaged) perturbed field is given by 

(B') = ( - Sro6 + ikzB^or‘>‘^r, B^q^^ -\- ik^B^r^^^^^, 

r or 


(23) 


It can be seen from Equation ([T4|) that includes a 


term of the form B.rQd^.^/dr, which woul d re sult in the in¬ 
clusion of a deriva tive o f in Equation (II) and thereby 
make the system (10 )-([l3| difficult to solve. The lead¬ 
ing coefficient in this case does not have a simple closed 
form that can be used to extract the resonances. There¬ 
fore, in the interest of maximizing our insight into the 
relevant resonances, we further simplify the treatment 
by letting Bpr be zero. Since the appearance of a radial 
field component is tied physically to the development of a 
vertically-sheared radial velocity field in the disk (in par¬ 
ticular, when Bz and dvorjdz are nonzero; see Equation 

t ), this approximation is consistent with our neglect of 
; VQr terms. 

With these two simplifications, the real part of the 
leading coefficient A 2 (r) of the resulting second-order dif¬ 
ferential equation has two obvious roots, which yield two 
resonance conditions. The first one is 


A-HZ. 

z' A 


m?a‘^ = 


>Az9rp^^p{kz-\-^)-\-dv^^c 




(24) 

Here vaz is the (rescaled) Alfven velocity associated with 
the z component of the field, = -0, and dv\^ = 

Ignoring order-h^ terms, the above 
equation indicates that one type of resonances occurs at 
the locations (one interior and one exterior to the planet’s 
radius) where the “corotating” perturbation frequency 
ma is equal to the frequency of a SMS wave propagating 
along Bz- As in Paper I, we refer to these as “Magnetic 
Resonances” (MRs). The second resonance condition is 


22 722 o 2 

TO cr =k^VAz- OVAy, 


1 - 




m 


3r2 


(25) 


^ Note that does not represent the midplane value of 
which is zero in our model. 
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If no t for the vertical gradient in the (p component, Equa¬ 
tion (251 would indicate that there is another type of res¬ 


onance at the locations where the perturbation frequency 
in the corotating frame matches that of a shear Alfven 
wave propagating in the z direction. As in Paper I, we 
refer to these as “Alfven Resonances” (ARs). Equations 
(24) and ([2^ indicate that the MRs and ARs for the field 


tio is generally > 1 near the midpla ne of a magnetically 
well-c oupled, wind-driving disk; e.g., [Konigl &: Salmeron 


configuration considered in this paper differ from those 
derived for the multicomponent field explored in Paper I 
(which is distinguished from the one discussed here by 
having a nonzero azimuthal component at the disk mid¬ 
plane) . We find that the MRs depend strongly on the ver¬ 
tical field and very little on the azimuthal field gradient 
(which only appears in conjunction with factors of order 
h^), whereas for the configuration considered In Paper I 
the MR locations are sensitive also to the azimuthal field 
component (see Equation (1.13)). Furthermore, in the 
limit oi kz 0 and no vertical field, the ARs are defined 
by m^cr^ = —dv\^{l — Thus, the ARs in this case 


only exist for values of m that satisfy 


u2 2 

h m 
3r2 


> 1, in stark 


contrast to the case of an azimuthal field with kz 0 
considered in Paper I (see Equations (1.11) and (1.12)), 
in which the ARs always exist. However, the ARs disap¬ 
pear, irrespective of the field geometry, in the strictly 2D 
limit (see Section 1.2.1). The 2D cases of dBjJdz ^ 0 do, 
however, exhibit weak MRs (see Equation (24)), which is 
consi stent with the r esults for a uniform azimuthal field 
(e.g., Terquem||2003 1, but they turn out to have a mini¬ 


mal effect on the net torque. 

The dependence of the resonance locations on the field 
configuration and on the wavenumber kz is illustrated in 
Figurej^ The MRs are shown by the black lines, and the 
ARs by the gray ones. For a pure vertical field, the ARs 
are always located outside the MRs (i.e., farther away 
from the corotation radius, where the disk’s angular ve¬ 
locity is equal to that of the planet). Both resonances 
have locations that depend on the azimuthal wavenum¬ 
ber TO, getting closer to the planet with increasing to. 
For a pure-R^ field, both of these resonances move fur¬ 
ther away from the planet as kz increases, diminishing 
the contributions of the resonances to the torque. For 
the B^o 7 ^ 0, Bzo = 0 case and kz ^ 0, the MRs are 
located symmetrically with respect to the corotation ra¬ 
dius, and their radial position does not vary with to. 
This is similar to the behav ior of the MR r esonances for 
the 2D, pure-R,^ case (e.g., Terquem 2003), except that 


in the case considered here tire MRs only appear away 
from the midplanej^ In this case, an additional AR ap¬ 
pears interior to the MR. As noted above, this AR only 
exists for wave numbers to above a critical value that de¬ 
pends on the disk thickness. For a thin disk with h 1, 
this critical value is large, which makes the torque con¬ 
tribution of these ARs negligible. 

Figure^ shows the outer (r > 1) wave propagation re¬ 
gions ana their relation to the locations of the turning 
points and resonances for the case of a pure vertical field 
and a finite value of kz. The magnetic field strength is 
parameterized by the value of j3z = I where (3 is 
half the thermal-to-magnetic pressure ratio. (This ra- 

® This can be seen in the numerical simulations presented in 
Figure 1.15. We capture the effect of these MRs even though our 
analysis is restricted to the midplane, where Bip is zero, because 
we use a vertically averaged value of 


2011|) Each of the resonances (MRs and ARs) is associ¬ 


ated with a pair of turning points that straddle it. We 
label the turning points according to the resonances they 
surround. From the planet’s location outward, there are 
four such turning points, Rm-) Rm-i-, Ra-, and Ra+, 
and then an additional one, Rl-i-) the modified effective 
outer Lindblad resonance. Two regions of wave propa¬ 
gation, from which torque is exerted on the planet, sur¬ 
round each resonance and are bounded by the turning 
points. Waves also propagate away from the outer Lind¬ 
blad turning point. The integrated torque exerted on the 
planet from the region r > 1 depends on the properties 
of these three regions as well as on the point-like con¬ 
tributions from each of the resonances. The net torque, 
in turn, is determined by the balance between the op¬ 
posing effects of the outer and inner disk regions, and 
by the contribution of the co-orbital region. Note that, 
as the value of the azimuthal mode number to increases, 
the resonances move closer to the planet but their asso¬ 
ciated regions of wave propagation become narrower. It 
is thus not obvious a priori which mode numbers dom¬ 
inate the torque. We address t his i ssue on the basis of 
o ur numerical results in Section 14.21 
Muto et al.j (|2008|) studied the pure-R^, kz 0 case 


using the shearing-sheet model and the WKB approxi¬ 
mation. Similarly to what we find, they identified three 
regions of wave propagation, corresponding (in order of 
increasing distance from the planet) to the SMS, Alfven, 
and FMS waves. Their Figures 1 and 2(c) (for which 
Pz = 0.9) are evidently most closely related to the re¬ 
sults we show in Figure Their inner evanescence re¬ 
gion, whose outer boundary (which they term LRU corre¬ 
sponds to our Rm- , exists only for k‘lv\^ > 3 Up PI In our 
work we also find that this region is present oruy when 
k1v\^ is sufficiently large, although the lower bound that 
we obtain with our more general treatment (3.04 D^) is 
slightl y different from the analytic limit. [Muto et al.j 
(20081 identified the inner boundary of the intermedi- 


ate evanescence region with the locus of the magnetic 
resonances rather than with that of the turning points 
Rm+ (which their derivation misses); we note in this 
connection that, even though we clearly distinguish be¬ 
tween these two loci, in practice they are very close to 
each other in our example for all values of to. The inner 
boundary of the intermediate evanescence region, given 
by the turning points Ra-j appears to correspond to 
what these authors call “vertical resonances.Their 
derivation, however, also misses the turning points Ra-i-, 
and thus they identify the inner boundary of the outer 
evanescence region with the locus of the Alfven reso¬ 
nances. In our example, the Alfven resonances are spa- 

® As was noted by [Muto et al.j (|2008j, this conditio n is the 
same as the linear stability condition against the MRI (e.g.,|Balbus| 

[2MT|. _ _ _ 

’^Note that [Muto et al.j ([2008^ refer to genuine resonances as 
well as to turning points as “resonances.” It is th us w orth reem¬ 
phasizing that the locations where K (Equation | |22^ ) vanishes, 
which include the classical effective Lindblad resonances in an un¬ 
magnetized disk, are actually turning points: They mark the onset 
of wave propagation in the disk, and, in contrast with actual res¬ 
onances, they do not correspond to a divergence in the perturbed 
density. 
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Fig. 1.— Locations of the magnetic (black lines) and Alfven (gray lines) resonances for two different field configurations, dBipJdz = 
0, Bz ^ 0 and dB^jdz yf 0, Bz =0 as a function of the azimuthal mode number m. The resonances for Bz yf 0 are shown for both 
kzh = 1.56 and kzh = 3.12 (the solid and dashed lines respectively), but the dB^jdz yf 0 case (dotted lin es! is dis played at only one 
vertical wavenumber since the resonance locations for this latter case are independent of kz (see Equations and | |25^ ). For the pure-i ?2 
case, the resonances move away from the planet as the vertical wavenumber increases, implying that the contribution of these modes to 
the net torque decreases with kz- The resonance locations for a field configuration with both Bz yf 0 and dBtpJdz yf 0 are very similar to 
those for the pure-B^ case, so this configuration is not shown separately. The Alfven resonances for the dB^ jdz ^ 0 case only exist for 
m > 55. 


tially well separated from the turning points Ra+ ex¬ 
cept at low values of m, where the WKB approximation 
(which requires m kzv) applies. This comparison il¬ 
lustrates the limitations of the WKB approximation in 
providing an accurate description of the disk-planet in¬ 
teraction. The outer boundary of the outer evanescence 
region is, however, correctly identified by this approxi¬ 
mation with the locus of the modified effective Lindblad 
resonances. 

In the limit = 0 of the pure-i ?2 case, both resonances 
vanish and there is only one set of turning points — the 
effective Lindblad resonances modified by t he presence 
of the field, as described by Equation (37) of|Muto et al. 
(2008). As noted above, our formalism does not involve 
the shearing-sheet approximation that was employed by 
these authors, and we are therefore able to evaluate the 
difference between the torques exerted by the inner and 
outer disk regions and hence the net torque acting on the 
planet. Figurej^shows the distances of the turning points 
on the two sides of the planet and their dependence on 
the field strength. It is seen that the basic structure is 
similar to that in a purely HD disk in that the outer effec¬ 
tive Lindblad resonances, which exert a negative torque 
on the planet, lie closer to the planet than the inner turn¬ 
ing points and can therefore be expected to dominate the 
interaction. This is true for all azimuthal mode numbers, 
although the asymmetry is stronger for low values of m. 
Under these circumstances, the planet is predicted to mi¬ 
grate inward. As the field strength increases, the turning 
points move farther away from the planet, which suggests 


that the magnitude of the net torque, and hence the rate 
of inward migration, are reduced. This heuristic infer¬ 
ence is borne out by our numerical results (see Figure]^ 
below) as well as by the 2D numerical simulations pre¬ 
sented in Paper I (see Figure 1.3). We further verified by 
numerical simulations that the decrease in the net torque 
with increasing applies also in 3D. 

When both B^o and Bzo are nonzero, the coefficients 
Aq, Al, and A 2 , and hence K (Equation (22)) are com¬ 
plex. In this case, we expect wave propagation when 
^{K} > 0 and iSi{K} :$> |S|jA|| as well as when 
|9{K}| |3fi{Ar}| (see Appendix]^. Figurej^shows the 

predicted regions of wave propagaHon (at a given m) for 
this field configuration for both 2D and 3D modes. The 
kz ^ 0 modes behave similarly to the 3D pure-H^ case, 
which can be understood from the fact that, even at the 
disk surface (subscript ‘s’), the amplitude of B^ remains 
a small fraction of Bz i\B^s/Bzo\ — 0.13 in this exam¬ 
ple). However, when kz = 0, there is an interior wave 
propagation region that starts very close to the planet 
(corresponding to a region of large |3{Ar}|), which dif¬ 
fers from the behavior of a pure-i?^ configuration in 2D. 
Although Figure]^ only shows the wave propagation re¬ 
gions for a given value of m, we have verified that the 
area of the interior propagation region in both the 2D 
and 3D cases decreases with increasing m (much as it 
does in the 3D pure-H^ model shown in Figure]^. The 
presence of an inner wave-propagation region in 2D for 
this field configuration suggests that, in contrast with 
the pure-i ?2 disk, the torque would be enhanced (rather 
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than reduced) in this case in comparison with the HD 
limit (see Section 4.4). 


2 . 6 . Effect of Vertical Gradient of B^ 

As we already pointed out, a vertical gradient of the ra¬ 
dial field component cannot be self-consistently incorpo¬ 
rated into our ideal-MHD disk model. However, in view 
of the fact that the magnetic tension force (oc BzdBr/dz) 
could play an im portant dynamical role in real wind¬ 
driving disks (e.g., |Wardle fc Konigl||1993 ), it is of inter¬ 
est to at least look for clues to its pos sible effect on planet 
migration. As we noted in Section |2.5[ when B^q 7 ^ 0, 
one canno t ob tain a differential equation in the form of 
Equation (20). We therefore attempt to apply an iter¬ 


ative approach to the procedure that yielded this equa¬ 
tion. First, we ignore the Brod^ip/dr term in Equation 
( [TT| ) and solve for as a function of fr and its derivatives 
only. We then differentiate the resulting expression for 
and use it to evaluate the term we originally ignored 
so as to obtain an updated value for We solve the 
rest of the system of equations using this updated value. 
In general, the leading-order coefficient of Equation (20 1 
that is obtained in this way is rather complicated, and 
there is no simple form for its roots. However, one can 
obtain approximate resonance conditions by keeping the 
dominant terms in the limit B^q —> 0 and also dropping 
the smaller-order cx h terms. In this way we find the 
following approximate conditions: 


99 9,9 „ dft / dr 

m a - vXzk^ - 2vAzdvAr - 


= 0 


m cr {vaz + C ) - VazC k, + dvArivAz + C )- 


dvA, 


dvAz / 2 2 


dr 


{vaz + C ) — 0 , 


(26) 


(27) 


where dvAr is defined analogously to dvAp in Section|2.5 


In the absence of a radial field. Equation (26) reduces 


to the Alfven resonance condition and Equation (27) to 
the magnetic resonance condi tion for a pure-i?^ config¬ 
uration. However, Equation (26) indicates that, when 


dvAr 7 ^ Oj there is an inherent asymmetry between the 
inner and outer Alfven-like resonances. This equation 
is formally cubic in cr, and since dfl/dr is negative, it 
yields an additional resonance at r > 1. This extra reso¬ 
nance persists even when kz = 0. This suggests that the 
presence of a vertical gradient of Br could enhance the 
torque exerted by waves exterior to the planet’s orbit, 
which induces inward migration. Although we cannot 
check on the validity of this prediction in the present 
study, we hope to pursue this case further in a future 
investigation. 

3. NUMERICAL SCHEME 

To solve Equation (|^, we employ the same methodol¬ 


ogy us ed by KorycansEy & Pollack (1993) and Terquem 
(2003). To save computational tinie, we employ known 


recurrence formulas to evaluate the Laplace coefficients 
(Equation (1^) for a given m. The relations we use are 
outlined in t he appendix of|Korycansky fc Poll ack|(|I993|) 
and a lso in Brouwer & Clemencell 96 1|) and Hagihara 

(p^. 


The integration of Equation (20) is generally tricky due 


of the various resonances. To achieve the required ac¬ 
curacy while using a standard initial-value integration 
routine, it is necessary to adopt small steps and require 
high precision; here we apply a fifth-o rder Runge-Kutta 
method with adaptive stepsize control ( Press et al.|1992 ) 
and use 16 -bit precision. To deal w i th any poles at a = 0 , 
we fol low Korycansky & Pollack (1993) and Terquem 
(2003) and displace the pole troin the real axis by id 
where 5 is a real number with |(5| ^ 1, i.e., a ^ a — iS. 
This effectively also displaces the resonances from the 
real axis. We start the integration from the planet and 
integrate toward the inner and outer domain boundaries 
(specified below). We find that this methodology pro¬ 
duces a smoother run than a “monotonic” integration 
from the inner boundary to the outer one. 

At the boundaries of our integration domain, where 
^ c^, we neglect Ai in comparison to Aq, which 
is justified by the fact that, to leading order, Aq/Ai 
{ r^j— fff j In this limit, the homogeneous 
version of Equation (pOl becomes 


+ Aoir)fr = 0 , 


(28) 


where we retain only the leading-order term in each co¬ 
efficient. Applying the WKB appro xim ation by assum¬ 
ing that the solution of Equation (28) takes the form 
yields few = A 0 /A 2 . For the case of a field 
'z, Aq is given by 


with nonzero Bz and dB^/d. 


An = 


22 2 

m a — K 


klc^K^ 


-klc^ 


(29) 


to the steepness of the planet potential and the presence 


Although the terms that include kz are small in com¬ 
parison wi th t he other terms on the right-hand side of 
Equation ( |2^ , they can influence the results at small 
values of m. Note also that, in the limit B —>■ 0 and 
kz —>■ 0 , A 0 /A 2 reduces to the dispersion relation for 
acoustic density waves in a 2D HD disk. 

We use this WKB result as a boundary condition, 
d^r/dr = ik^^^r, at the inner and outer edges of the 
integration domain. We still need initial values for 
and dfridr to start the integration, so, similarly to the 
standar d shoo ting method and follo wing [Korycansky fc| 
Pollackl (|1993 ) and Terquem (2003), we start with ran¬ 
dom values and use an iterative approach to capture the 
true solution. For the sake of completeness, we provide a 
summary o f this procedure below; additio nal details can 
be found in Korycansky & Pollack (1993). 

Starting with random values implies that the solution 
obtained by the numerical scheme will in general con¬ 
tain a linear combination of solutions of the homoge¬ 
neous equation (subscript ‘hs’) and of the particular so¬ 
lution of the inhomogeneous equation (subscript ‘ps’), 
i.e., fr = fps + C'lChsi + C' 2 ^hs 2 - To attain conver¬ 
gence toward the true solution, we simultaneously in¬ 
tegrate the two linearly independent homogeneous so¬ 
lutions, ^hsi and ^hs 2 - After one full integration, we 
can use the boundary conditions to constrain the con¬ 
stants Cl and C 2 (specifically, our integrated solutions 
for ^ps, Chsi; Chs 2 and their derivatives yield a 2 x 2 system 
for Cl and C 2 when plugged into d^r/dr = iki^rfr at the 
disk boundaries). After solving for these constants, we 
perform one more integration, this time initializing ^ps 
as 5 ps,oid + C'lChsi + C' 2 ^hs 2 — this allows ^ps to converge 
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Fig. 2.— Wave propagation (clear) and evanescence (shaded) regions outside the planet’s orbit for a purely vertical field configuration 
with 0z = 0.8 as a func tion of the azimuthal mode number m. The vertical wavenumber is fixed at kzh = 1.56. The magnetic resonance 
given by Equation | |24[ | is labeled ’MR’, whereas the Alfven resonance given by Equation | |25| | is labeled ’AR’. The turning points are 
labeled Ra+-, ^A- > ^M+? Rm— (see text for details). 



Fig. 3.— Wave propagation (clear) and evanescence (shaded) regions on both sides of the planet’s orbit for a purely vertical field in the 
2D limit. The three panels correspond to different magnetic field strengths, parameterized by I3z Note the asymmetry between 

the inner and outer turning points (dashed lines) in each case. 
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Fig. 4. — Distribution of wave propagation (clear) and evanes¬ 
cence (shaded) regions in a disk with a Bz ^ 0, dB^jdz ^ 0 field 
configuration for kzh = 0 (top) and kzh ^ 0 (bottom). The 2D 
case corresponds to the same parameters as the left panel of Fig¬ 
ure]^ although in this schematic we neglect the contribution of the 
order-h^ terms. The bottom panel corresponds to the kzh = 3.12, 
m = 20 case shown in Table [T] The dashed lines indicate the 
locations of the Alfven and magnetic resonances. 


toward ^j.. As in Terquem (20031, we find that good 
convergence is attained already atter two iterations. 

In contrast with the behavior of the 2D, purely az¬ 
imuthal configuration considered by Terquem (2003), in 
which, for example, the location of the maghetic reso¬ 
nance was independent of to, for the field configurations 
that we consider the locations of the resonances vary with 
the azimuthal mode number to as well as with the vertical 
wavenumber kz- In our case, large vertical wavenumbers 
and small azimuthal mode numbers correspond to reso¬ 
nances (and, more generally, wave propagation regions) 
that lie far from the planet. Calculating the torque can 
then be computationally expensive, as it requires integra¬ 
tion over a large region to capture the full effect of the 
resonances. To overcome this challenge, we limit our¬ 
selves to low values of kzh (although we also calculate 
the torque for a few higher values of this product) and 
interpolate between some of the larger values of to. For 
very high values of to, we find it necessary t o re duce 
the size of the integr ation domain (see Sectio n 4.4 ). As 
was also the case in Korycansky & Pollack (1993), we 
find that for certain disk parameters we cannot calculate 
the contribution of some low azimuthal mode numbers 


(the differential equation becomes stiff and the needed 
step size to integrate is excessively small). In these cases 
we extrapolate to estimate the torque contribution from 
these low to. This, in turn, introduces an error into the 
summed (over to) torque value, but we find that the qual¬ 
itative behavior of the sol ution is not affected by this 
pro cedure. Note tha t, as in|Korycansky & Pollack|( 19931 
and|Terquem| (|2003|) , we only attempt to cafcufate modes 
with TO > 2. For to = 0, there is no effect on the torque 
since oc to (see Equation (30)), whereas for to = 1, 
the resul ts involve a con tribution from the disk’s outer 
edge (see|Shu et al.|jl990|, which we eliminate by enforc¬ 
ing an outgoing wave co nditi on at the radial boundaries. 
Solutions to Equation (20) are highly oscillatory, lead¬ 


ing to strong cancellations of the torque contributions 
from adjacent zones, particularly in the outer regions 
of the disk. For this reason, the numerical results can 
depend on the size of the integration domain. In our 
treatment, we use Hn = 0.3 Tp and rout = l-7rp as the 
integration boundaries for to > 6. This domain is slightly 
larg er than the one u sed by Korycansky & Pollack (1993) 
and Terquem (2003) since we need to ensure that the in¬ 
tegration region is sufficiently large also for the kzh ^ 0 
cases, in which resonances and turning points lie farther 
away from the planet. For to < 6 we adopt an even 
larger domain, r-m = 0.2 rp and rout = 5.0 rp. We assume 
that the torque value for a given to has effectively con¬ 
verged if the value does not differ by more than ~ 1% 
from the torque obtained using a slightly larger range . 
Foll owing the pract ice in Korycansky & Pollack (1993) 
and Terquem (2003), we set the values of the parameters 
S and e to be 1 X 10“® and 1 x 10“^, respectively; how¬ 
ever, our results are not sensitive to the specific choices 
of these values provided that they remain small (either 
one can be a factor of ^ 10 larger before any changes are 
detected in the results). 

After solving for ^r, we obtain W' from Equations (10) 


- (13) and proceed to calculate the torque (at a given to) 
per unit area (A) that is exerted by the planet using 


dT 

u±m 

~Wa 


- \ iWY^oip' 


and 


Tm = 27r 


dTrr. 

dA 


rdr 


(30) 


(31) 


(see Korycansky & Pollack 1993). To get the cumulative 
torque, we sum over to (which iii practice we do through 
an integral, assuming that a given value of Tm charac¬ 
terizes a small interval of to values). To obtain the total 
torque in the 3D cases, we approximate the additional 
sum over the vertical wave numbers by adding up the 
contributions of the specific kz values that we evaluate. 

In this paper we use the convention, commonly 
adopted in linear analyses of the type we ca rry out (e.g., 
Korycansky & Pollack|l993 Terquem 2003), wherein one 
calculates the torque exerted by the planet on the disk. 
In this case, a positive torque implies inward migration 
whereas a negative torque indicates outward migration. 
This sign convention is the opposite of that commonly 
used in numerical studies, which concentrate on the effect 
that the disk has on the planet. The latter convention is 
also adopted in Paper I, where the focus is on numerical 
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simulations. We caution the reader to be mindful of this 
difference when comparing the results of the two papers. 

4. RESULTS 

4.1. Hydrodynamic Limit 

We ran ou r scheme with = 0 to che ck against the 
HD results of Korycansky & Pollack (19931. In this limit, 
waves propagate away from the ettective Lindblad reso¬ 
nances and there is a genuine resonance at the corotation 
radius. The total torque exerted by the planet is posi¬ 
tive, representing the standard inward-migration result. 
We found t hat ou r results reproduce those of |Korycansky| 
'& Pollack| (|1993|) well. In particular, we obtained a good 
match to their plot (in Figure 2) of the real and imaginary 
parts of W{r) for m = 10, and our calculated cumula¬ 
tive dimensionless torque of ^ 1300 is in good agree¬ 
ment with the value reported in their Figure 14. Our 
results similarly match the HD calculations presented in 
Terquem (2003), and we further validated our numerical 
scheme by reproducing the azimuthal field results given 
in that paper (in particular, the plot of Tm vs. m for the 
uniform-field case, shown in the top panel of Figure 6 of 
that paper) p] 

We also c onsidered the 3D torque for this case. |Tanaka| 
et al.|(|2002|) originally found that the 3D torque is weaker 
than the 2D torque by a factor of > 2. For the kzh values 
used in this paper, we found the ratio to be closer to ~ 1. 
We note, howe ver, that an exact co rrespondence is not 
expected since |Tanaka et al.| (|2002|) used different disk 
parameters and included vertical stratification. Further¬ 
more, the value of the 3D torque in our analysis also de¬ 
pends on the range of k^h that we use (which in practice 
corresponds to the choice of h, the effective disk half¬ 
thickness). Our chosen minimum value of kzh (= 1.56) 
is based on applying our adopted value of fiz (= 0.8) to 
the condition for the existence of the inner evanescence 
region (or, equivalently, the condition for the existence 
of the i nner most turning point) in the pure-Hj, case (see 
Section 2.5). Regarding the upper limit of the adopted 


a factor ~ 2 sm aller than the HD value. As was alre ady 


range, it is m practice limited by the computational diffi¬ 
culties encountered for large values of kzh (See Section!^. 
We were, however, able to explore the 3D torque in tnis 
case for sufficiently large values of this product (repre¬ 
sented by ever thinner disks) to confirm that it is indeed 
weaker t han the 2D HD torqu e, in qualitative agreement 
with the Tanaka et al. (2002) result. 


4.2. Bz ^ 0, dB^/dz = 0 Field Configuration 

Although |Muto et al.| (|2008D investigated the case of 
a pure-i ?2 held, they formulated the problem in the 
shearing-sheet approximation and considered the 3D 
regime only for disks that are strongly magnetically dom¬ 
inated {fiz < 10“^) and thus do not correspond to wind¬ 
driving systems. We therefore take a fresh look at this 
case. 

Figure shows the integrated torque as a function 
of m for a uniform disk in the 2D limit and for two finite 
values of kz- We also plot the HD result for compari¬ 
son. The 2D integrated torque is invariably lower than 
in the unmagnetized case, resulting in a cumulative di¬ 
mensionless torque for the chosen disk parameters that is 

The 2D, pure-B ,,3 case was calculated using the formulation 
presented in Section 1.2. 


pointed out by Muto et al. (2008) (see also Section 2.5 
and Paper I), this reduction can be attributed to the 
outward shift of the effective Lindblad resonances in a 
magnetized disk, which increases with the field strength 
(see Figure]^. The behavior of the 3D modes is more 
complicated. For kzh = 1.56, is positive and gener¬ 
ally larger than the corresponding value for kz = 0. This 
is due to the appearance in this case of the magnetic and 
Alfven resonances and of their associated regions of wave 
propagation. However, for the larger vertical wavenum¬ 
ber considered in this figure {kzh = 3.12), the net torque 
for each value of m is negative. This is a consequence of 
the fact — revealed by a close examination of Figure]^ 
— that the inner MR lies slightly closer to the planet 
than the outer MR. This asymmetry becomes more pro¬ 
nounced, with a corresponding increase in its effect, as 
kz increases. We thus expect the integrated torque to 
remain negative as kzh grows further; however, its con¬ 
tribution to the total torque would be smaller than that 
of the kzh = 3.12 mode because the resonances and asso¬ 
ciated wave-propagation regions move outward with in¬ 
creasing kz ■ We estimate the total torque in this case by 
adding up the net torques for the three values of kz shown 
in Figure]^ we find a value that is still smaller (^ 0.7) 
than the 2D HD case shown in the same plot. This ball¬ 
park estimate is in good agreement with the total torque 
obtained in the numerical simulations presented in Pa¬ 
per I (see Figure 1.10). In Section 3.4.1 of that work we 
point out that, for the adopted model parameters, our 
semianalytic estimates yield comparable values for the 
3D and the 2D HD torques. Hence, our conclusion that 
the torque in a disk with a pure-Hz held is roughly half 
that of an unmagnetized disk applies in both 2D and 3D. 

The breakdown of the contributions to the integrated 
torque from different regions in the disk for kzh = 3.12, 
fiz = 0.8, and two values of m (20 and 100) is shown 
in Table For m = 20, we also show results for 
a field that is twice as strong {fiz = 0.4). The re¬ 
gions are chosen as follows (see Figure for the mean¬ 
ing of the different labels): from the inner/outer edge 
of the disk to the inner/outer Rl+ {T^i and T^^o, re¬ 
spectively); from the inner/outer Rl+ to the inner/outer 
Ra- {TARi and Taro, respectively), and from the in¬ 
ner/outer Ra- to the planet’s location {TMRi and Tmro, 
respectively). These sectors encompass the propagation 
regions of FMS, Alfven, and SMS waves, respectively. 
Our ability to distinguish between the contributions of 
the regions on the inner and outer sides of the planet 
helps us gain insights that could not be obtained with a 
shearing-sheet model. 

The comparison between the m = 20 and m = 100 
results for fiz = 0.8 reveals that the contribution of the 
magnetic resonance region dominates, and its magnitude 
increases with m, on each of the two (inner and outer) 
sides of the planet. However, the contributions from the 
two sides are closer in magnitude, resulting in a smaller 
net torque, for the higher-m case. On the other hand, 
the contribution from the Alfven resonance region (from 
either side of the planet as well as the net one) decreases 
with m. These two trends combine to make the inte¬ 
grated torque much higher (with comparable contribu¬ 
tions from the magnetic and Alfven resonance regions) 
for m = 20. This behavior also explains why peaks 
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at a comparatively low value of m ( see Figure]^. 

The torque arising from the regions around a resonance 
is a combination of the point-like contribution from the 
resonance location and that from the associated wave- 
propag ation zone . In th e case of the MR we find, simi¬ 
larly to Terquem (2003) for the pure-R,^ field, that the 


point-like contribution to Tm typically has a smaller mag¬ 
nitude than the contribution from the surrounding re¬ 
gion, and that these two contributions can have opposite 
signs. (For the most part, the direction of the point-like 
torque oscillates at low values of m and converges toward 
a set sign at high m.) The same behavior holds true for 
the AR, but in general the point-like contribution from 
the AR is much smaller in magnitude in comparison with 
both the contribution of the surrounding AR region and 
the point-like MR torque. 

Table also shows the breakdown of the integrated 
torque for a field that is twice as strong, i.e. /3^ = 0.4, 
at TO = 20. Equations (24) and (25) indicate that, when 


the field becomes stronger, the magnetic and Alfven reso¬ 
nances move further away from the planet, which has the 
effect of decreasing the magnitude of the planetary po¬ 
tential at their locations. We find, however, that while 
the strength of the disk response (as measured by the 
magnitude of 3{1U'}) indeed decreases in the AR region 
(as well as in the FMS wave propagation region beyond 
the outer Lindblad turning point) when B^o is increased, 
other factors conspire in this case to increase the enthalpy 
amplitude in the MR region. We have verified this trend 
with numerical simulations, which show a more defined 
wake from SMS waves but a weaker Lindblad wake when 
j3z is decreased. The decrease in torque from the Alfven 
and FMS wave propagation regions even as the contri¬ 
bution from the MR region goes up results in a slight 
reduction in the total torque when 15^ is decreased from 
0.8 to 0.4. Since the torque at a given value of to is not 
always representative of the cumulative torque (typically 
on account of an ancillary effect such as the variation in 
the sign of the point-like contribution described in the 
preceding paragraph), we calculated the torque for the 
(3z = 0.4 case over the same range of to and k^h that was 
employed in Figure We found that the total torque 
calculated in this way was again slightly smaller for the 
stronger field. We also carried out numerical simulations 
for the two values of /3z considered in Table and con¬ 
firmed that the torque is reduced for the higher value of 
Bz- 

The numerical estimates given in the caption to Fig¬ 
ure indicate that, for /3z = 0.8, the total torque on the 
planet is roughly the same in 2D and 3D. If this infer¬ 
ence is correct, the implied behavior would contrast with 
that of HD disks, in which the 2D torqu e dominates (e.g., 
Tanaka et al.|2002 ). Muto et al. (2008) inferred that the 
3U torque exerted on a given side of the planet exceeds 
the corresponding 2D torque in the limit /3z <C 1 (for 
which the MR contribution dominates), but this regime 
is not relevant to the formation region of giant planets. 
Furthermore, their analysis was carried out in the con¬ 
text of the shearing-sheet approximation and therefore 
could not anticipate our findings (for /3z < 1) that the 
3D torque decreases only slightly with increasing field 
strength for given values of kzh and to because of the 
opposing effects of the two sides of the planet, and that, 
in addition, the kz ^ 0 contributions to the total torque 


largely cancel each other outj^ 

4.3. dB^jdz ^ 0, Bz = 0 Field Configuration 

We consider this field configuration both here and in 
Paper I even though it is not realistic in an attempt to 
isolate the effect of the dB^p/dz term on planet migration 
in wind-driving disks. Figure|^shows the net torque Tjn 
vs. TO for this configuration in the 2D limit and for the 
two vertical wavenumbers considered in Figure As 
in the pure-R^ case shown in the latter figure, we were 
unable to derive results for a few 1ow-to modes in the 
kzh — 3.12 case, but their effect is likely not crucial in 
this case either. We find that the contributions from 
the two fcz ^ 0 integrals that we present roughly cancel 
each other out as they did in the pure-Rz case, and we 
thus approximate the total torque by the contribution of 
the kz = 0 mode. In this way we infer a value for the 
2D torque that is nearly equal to that in the HD limit, 
which agrees with the result of the numerical simulations 
presented in Paper I (see Fi gure 1.10). 

As we discussed in Section [231 MRs appear above and 
below the midplane in the 3D solutions for this case, and 
the evidence for them — seen in the behavior of the per¬ 
turbed enthalpy — persists into the 2D limit. In this 
limit, the MRs can be regarded as the “smeared out” 
(over the vertical extent of th e disk) version of the 2D, 
pure-R,p resonances found by Terquem (2003). These 
“diluted” MRs are, however, quite weak, and their con¬ 
tribution to the torque is very small. In addition, there 
are no ARs for this field configuration in the 2D limit. 
This is because the terms in the coefficients of Equation 
(20) that would have given rise to these resonances can¬ 
cel out in the limit kz, Bz —>■ 0{^ These facts explain 
our finding that the magnitude oithe 2D torque for this 
case is close to its HD value. 

We also already pointed out (in connection with Equa¬ 
tion (25)) that, in 3D, ARs only exist in this case for 


TO^ > 3r^//i^. When this condition is satisfied, the ARs 
appear very close to the planet, even more so than the 
MRs (see Figure[^. The larger the value of h, the smaller 
the critical value of to, and since for low values of to the 

f ' inet’s jwtential couples well to the disk (see Equations 
and Q), the stronger could be the influence of the 
met onrhe disk. In the example presented in Figure [^ 
we chose the disk half-thickness to be of the order of what 
would have been the density scale height (c/U) had we 
explicitly included the vertical stellar gravity. This was 
sufficient to guarantee that, despite the fact that we as¬ 
sumed a uniform-density disk, the effect of the ARs on 
the integrated torque did not become artificially large. 

4.4. Rz 0, dBpjdz ^ 0 Field Configuration 

Having considered the effects of the Bz 0 and 
dBp/dz ^ 0 field configurations separately in the previ¬ 
ous two subsections, we now investigate the outcome of 


W e do, however, find evidence for the trend that |Muto et al^ 
l |2008[ l identified — that, as decreases, the evanescence regions 
grow and the point-like contribution of the MRs goes up — also in 
the higher-regime that is of interest to us (see Table 

This is not apparent from the results given in Section!^ where, 
to simplify the presentation, we do not write down the expressions 
for the coefficients of Equation ( |20[ |. However, in the related case 
of a pme-Bip field considered in t^aper I, we show explicitly that 
the ARs vanish in the 2D limit (see Equation 1.15). 
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Fig. 5.— Spatially integrated torque exerted by the planet on a thin {h = c/{rpQp) = 0.03), uniform disk that is threaded by a purely 
vertical field (of strength I3z = 0.8) as a function of the azimuthal mode number m. Results for the 2D limit (kz = 0) and for two nonzero 
vertical wavenumbers are shown, with the 2D HD case included for reference. Because of the stiffness of the governing differential equation 
at low values of m, we were only able to obtain numerical results for m > 8 for kzh = 3.12, m > 4 for kzh = 1.56, and m > 2 for the 
kzh = 0 case. Extrapolating where necessary, the cumulative dimensionless torques for kzh = 0, 1.56, and 3.12 are ~ 700, 1400, and —1100, 
respectively. Neglecting the contribution from higher kz modes, we thus estimate the total dimensionless torque to be ~ 1000, a fraction 
~ 0.7 of the 2D HD value 1300). 


combining them together. As we discussed in connection 
with Equation the simultaneous presence of these 
two terms in the angular momentum conservation equa¬ 
tion gives rise to vertical angular momentum transport 
in a magnetically threaded disk, and they are thus a key 
ingredient of a wind-driving disk model. We parametrize 
this configuration throu^the values of dB^/dz (which 
equals B^q] see footnote^ and B^o (or, equivalently, of 
Pz and Pip). The ratio oirhe azimuthal field amplitude 
at the disk surface to the vertical field component can 
be expressed as \Bips/Bzo\ = h^/pz/Pip, and is equal to 
0.13 for our chosen parameters. This value is consistent 
with models of wind-driving disks that are magnetically 
active at the midplane, for which this ratio is t ypically 
inferred to be ^ 1 (e.g., Wardle fc K5nigl||l993 ). 

Not surprisingly, the behavior of the “combined” field 
configuration differs from that of its separate compo¬ 
nents. We noted in Section [231 that one distinctive fea¬ 
ture of adding the Bp gradient term to the uniform Bz 
term in the 2D limit is the emergence of a wave propaga¬ 


tion region in the immediate vicinity of the planet (see 
top panel of Figure |4). We illustrate the detailed struc¬ 
ture of this region inFigurej^ which plots the radial de¬ 
pendence of the perturbed enthalpy W' for two different 
values of Pp (while keeping the magnitude of Pz as well 
as the value of m fixed). It is apparent from the shape of 
the imaginary component o f W ' (which enters into the 
torque expression. Equation (301) that these solutions ex¬ 
hibit wave-like behavior near the planet, and that this be¬ 
havior becomes more pronounced as \BpQ/Bzo \ increases. 
This region can be expected to make a strong contribu¬ 
tion to the torque, and we indeed find (see Figurethat 
the integrated torque in this case continues to grow with 
decreasing m, and that, correspondingly, the cumulative 
torque is significantly larger than the value for either of 
the 2D cases shown in Figures and Another notewor¬ 
thy finding is that the contributions to the torque from 
this region have the same (positive) sign on both sides 
of the planet. This behavior, which was not previously 
encountered in planet-disk interactions, is also exhibited 
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by the kz ^ 0 modes for this field configuration (see next 
paragraph). 

As the bottom panel of Figure [^illustrates, the distri¬ 
bution of the wave propagation regions for the ^ 0 
modes of this field configuration is qualitatively similar 
to that of a pure-Uj, field. In fact, given the smallness of 
the factor h?dv\^ = « 0.02 in Equations 

(24) and (25), the locations of the MRs and of the ARs 
are essentially the same in these two cases. However, the 
disk response is very different. This is demonstrated by 
the last two lines in Table which provide the break¬ 
down of the contributions to the integrated torque in this 
case for the same two values of m as in the first two lines 
of this table, which list the corresponding con trib utions 
for the pure-Hz field configuration (see Section 4.21. It is 


seen that, whereas the contribution of the MR regions in¬ 
creases with m and that of the AR regions decreases with 
m in both cases, other major aspects of the disk behav¬ 
ior have changed. In particular, while the regions interior 
and exterior to the planet contribute with opposite signs 
in the pure-Rz case (as they also do in an HD disk), in 
this case the contributions from the MR and AR regions 
have the same sign on both sides of the planet for suffi¬ 
ciently high values of m {> 20 in this case): Those from 
the MRs are positive everywhere, whereas those from the 
ARs change from negative to positive as m increases (and 
we verified that in each case the point-like contribution 
has the same sign as that from the surroundings). Thus, 
in contrast with the behavior of the pure-Rz held, where 
the contributions from the two sides of the planet offset 
each other and lead to comparatively small integrated 
torques, in this case the growth of the MR torque with 
increasing m, and the fact that the contributions from 
the two sides of the planet add up, can result in very 
large integrated torques. 

Calculating the cumulative torque for the fcz yf 0 
modes in this case is complicated by the fact that the val¬ 
ues of Tm are still growing even as m approaches ~ 100 
(see Figure]^. Integrating dT^jdA at high values of m 
is computationally challenging because of the high oscil¬ 
lation frequency, which necessitates taking ever smaller 
step sizes. We have overcome this difficulty by reducing 
the integration-domain boundaries for high values of m, 
making however sure that, with each such reduction, the 
torque already calculated for a given value of m did not 
change. The inset in Figure shows the result of this 
calculation up to m = 1500 for the lowest vertical mode 
{kzh = 1.56). Since, for high values of m, the contribu¬ 
tion from the regions farthest away from the planet (the 
“Lindblad” and AR regions; Ti^i/o and TARi/o in Tablej^ 
are negligible, the value of the integral is not meaning¬ 
fully affected by shrinking the integration domain. At 
these very high values of m, the contribution from the 
MRs — which is added with the same sign from both 
sides of the planet - dominates. However, as the mode 
number continues to increase, the physical width of the 
resonance becomes progressively smaller (corresponding 
to a decreasing width of the wave propagation region, 
akin to the high-m behavior of the pure-Rz field seen 
in Figure and eventually (for m > 700) the point¬ 
like contribution becomes negligible and Tm converges to 
zero. The behavior of the other vertical modes is similar 
and they evidently also contribute a net positive torque 


(the case kzh = 3.12 is shown in Figure]^ for m < 100); 
however, since the contribution of the kzn~= 1.56 mode is 
by far the dominant one, we have not pursued the higher 
kz modes to full convergence. 

5. DISCUSSION 

The main goal of the present analysis has been to gain 
insights into Type I planet migration in wind-driving, 
magnetized disks. A key property of such disks is verti¬ 
cal angular momentum transport, which is brought about 
by the magnetic torque term oc BzdB^/dz. Given that 
the equilibrium structures of such systems cannot be 
studied self-consistently in the context of ideal MHD, 
we considered separately the behavior of pure-Rz and 
Y>VLTe-dB^pldz disks, and then attempted a linear pertur¬ 
bation analysis also for the combined Rz -I- dB^/dz field 
configuration. 

In the case of a uniform vertical field component, we 
determined that the torque is reduced in comparison with 
the pure-HD torque in both 2D and 3D. This implies that 
planet migration would still be inward, but at a lower (by 
a factor of ~ 2) rate than in a disk not threaded by a 
large-s cale field. For th e cas e of a uniform azimuth al field 
in 2D, |Terquem| (|2003|) and |Fromang et al.j (|2005D found 
that the torque is actually enhanced over tne HD torque 
(resulting in faster inward migration). However, these 
authors also demonstrated that if \B^\ decreases suffi¬ 
ciently rapidly with radius then inward migration can be 
strongly suppressed for {B^{r) oc r~^) and even reversed 
(for B^{r) oc r~^). The physical reason for this behav¬ 
ior is that an outward decrease in the field amplitude re¬ 
duces the contribution of the outer MR region (where the 
torque exerted by the planet is positive) relative to that 
of the inner MR region (which has the opposite effect). 
An outward-decre asing field amplitude (corresponding to 
Qz < 0; see Section jOf can be expected to arise naturally 
in the case of an open magnetic field that is dragged in¬ 
ward by the disk accretion flow, and it is thus of interest 
to inquire whether a similar reversal in the direction of 
planet migration is possible also for the astrophysically 
more relevant case of a vertical midplane field. We now 
show that the answer to this question is no. 

In the 2D limit, a stronge r ve rtical field leads to a 
decreased torque (see Section 4.21, in contrast with the 
pure Bcp case, where the converse is true. Thus, if Rz 
decreases with r, the torque from the inner region of the 
disk (r < 1) will be even weaker relative to the torque 
from the outer region (r > 1) than in the uniform-Rz 
case, resulting in a more positive overall torque and hence 
in faster inward migration. An added effect arises from 
the fact that a negative value of dBz/dr corresponds 
to an outward-decreasing magnetic pressure that acts to 
counter the stellar gravity, thereby shifting the corotation 
radius inward (see Equation (1.17)). This, in turn, places 
the outer magnetic resonances and associated turning 
points closer to the planet, resulting in a further enhance¬ 
ment of inward migration. In the 3D case, our numer¬ 
ical simulations and semianalytic results indicate that 
the total torque similarly decreases with increasing field 
strength for the range of /3z values we consider, and that 
the magnitude of the torque is comparable to that in 2D 
and is again lower than the HD torque. This leads us to 
infer that it is unlikely that a realistic radial gradient in 
Rz could give rise to outward migration. 
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Fig. 6.— Same as Figure]^ but for a dBi^/dz ^ 0, Bz = 0 field configuration. In this example, = 0.04. For the kzh = 3.12 case 
we were able to obtain numerical results only for m > 8, and, for kzh = 1.56, only m > 2. Extrapolating where necessary, we find the 
cumulative dimensionless torque for the kzh = 0, 1.56, and 3.12 cases to be ~ 1200, 1000, and -900, respectively. Thus neglecting the 
contribution from higher kz modes, we estimate that the total dimensionless torque is similar to the value for the HD case (~ 1300). 




Fig. 7.— Spatial behavior of the real (dashed line) and imaginary (solid line) parts of the perturbed enthalpy for the m = 10 mode, 
in the 2D limit of the Bz ^ 0, dB^pjdz ^ 0 field configuration. Results are shown for a fixed vertical field amplitude {^z = 0.8) but two 
different azimuthal field parameters: {Sip = 0.04 (left) and /3 = 0.001 (right, corresponding to a \/40 times higher value of \dBip/dz\). The 
small sharp peaks exhibited by the dashed curve near r = 1 correspond to the order-fi^ terms in the expression (Equation (|24[ |) for the 
MRs and have negligible effect on the torque. The normalized integrated torque Tio is ~ 100 for the left panel and ~ 300 for the right one, 
indicating that the magnitude of the torque increases with the strength of the field gradient. 
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Fig. 8.— Same as Figure]^ but for a dB^/dz ^ 0, Bz ^ 0 field configuration. In this example, ^z = 0.8 and 13^ = 0.04. We estimate 
a torque for the 2D kz = 0 mode of ~ 3400 in dimensionless units, nearly 3 times larger than the torque in the HD case (~ 1300). The 
kz ^ 0 modes converge to zero much more slowly with increasing azimuthal mode number m. The inset shows the result for the (dominant) 
kzh = 1.56 mode, which yields an integrated dimensionless torque of ~ 230, 000, a factor of ~ 200 larger than the HD torque. 


TABLE 1 

Net torques for different field configurations in 3D (kzh = 3.12) 


B 

/3z 


m 

TLi 

Tari 

TMRi 

TmRo 

Taro 

Tlo 

^integrated 

B, 

0.8 

-- 

20 

-0.5 

-16.5 

-40.7 

32.1 

11.4 

0.2 

-14.0 

Sz 

0.8 

— 

100 

-0.1 

-2.3 

-65.2 

63.6 

2.2 

0.1 

-1.7 

Sz 

0.4 

— 

20 

> -0.01 

-2.2 

-54.5 

43.0 

0.1 

< 0.01 

-13.6 

B:,+dB^ldz 

0.8 

0.04 

20 

-0.5 

-20.4 

15.3 

15.3 

-7.0 

0.1 

2.8 

B^+dB^/dz 

0.8 

0.04 

100 

-0.2 

0.2 

60.9 

59.7 

0.8 

0.1 

121.5 
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In the case of a pure dB^jdz field configuration, we 
found that the effect of the magnetic field on the planet- 
disk interaction is minimal. In fact, the behavior of such 
a disk is very similar to that of an HD disk, although it 
also exhibits weak magnetic (in 2D and 3D) and Alfven 
(in 3D) resonances away from the midplane. However, 
when both a vertical field and a growing azimuthal field 
component are present, we discovered that the torque 
exerted by the planet on the disk is greatly enhanced 
compared to the HD case. In the 2D limit, ou r tu rning- 
point analysis for this configuration (Section 2.5) indi¬ 
cates that wave propagation occurs in the vicinity of the 
planet. The torque on this close-in region dominates the 
interaction and leads to inward migration that is faster 
by a factor of ^ 3 than in an HD disk. The 7 ^ 0 
modes also contribute to inward migration, but the nu¬ 
merical evaluation of their cumulative torques is compli¬ 
cated by the fact that the integrated torque for each 
of these modes converges to zero only at very high val¬ 
ues of m (> 10^). We pursued a full calculation only 
for the lowest (but, by far, the dominant) vertical mode 
{kzh = 1.56), which yielded a cumulative torque that is 
a factor of ^ 200 larger than the HD value. In this case 
the torque is dominated by the contribution of the MRs, 
which move closer to the planet as m is increased. The 
origin of the comparatively strong torque in both the 
2D and 3D regimes is the fact that, unlike the situation 
in the HD and pure-i?^ cases, where the torque contri¬ 
butions from the resonance regions on the two sides of 
the planet add with opposite signs and partially cancel 
each other out, for this field configuration they combine 
with the same sign. This, in turn, suggests that a ra¬ 
dial gradient in the field strength would only affect the 
speed of the planet’s radial drift in this case, but not 
its (inward) direction. In Appendix [B| we show that the 
behavior of the disk for this field connguration is tied to 
the term oc B^dB^/dz in the angular momentum conser¬ 
vation equation. Thus, in this case the planet evidently 
taps directly into the vertical angular momentum trans¬ 
port channel that enables radial accretion (i.e., vor < 0 ) 
in the equilibrium disk model, and this becomes the main 
mechanism through which the planet loses angular mo¬ 
mentum and moves (very rapidly) inward. 


6. CONCLUSION 

This paper and its companion (Paper I) explore the 
effects of a large-scale, ordered magnetic field in proto¬ 
planetary disks on Type I planet migration. A large- 
scale, open field can be expected to be present in the 
planet formation regions of such disks on account of the 
advection by the disk accretion flow of the interstellar 
field lines that thread the parent molecular cloud core. 
The inward dragging of the field lines would generate a 
radial field component, and the latter would, in turn, be 
sheared by the differentially rotating disk material and 
give rise to an azimuthal field component. This would 
naturally give rise to a magnetic torque (oc B^dB^/dz) 
and hence to vertical angular momentum transport, pos¬ 
sibly in conjunction with the development of a centrifu- 
gally driven disk wind. Our primary motivation in this 
work is to determine the behavior of planets that reside in 
disk regions of this type. This task is challenging because 
quasi-steady equilibrium models of disk regions that are 
subject to vigorous radial advection and azimuthal shear¬ 


ing require the application of nonideal MHD, which is a 
complex undertaking. In preparation for a full treatment 
along these lines, we set out to investigate simpler model 
problems that can be treated using ideal MHD. In Pa¬ 
per I we examine a range of problems of this kind through 
a combination of numerical simulations and a linear per¬ 
turbation analysis. This paper complements that work 
by focusing on the semianalytic approach and applying 
it to the Bz + dB^p/dz field configuration, which, despite 
being at the heart of the wind-driving disk model, cannot 
be simulated using an ideal-MHD numerical code. 

While the calculations undertaken in this paper employ 
a number of simplifications, we present a general formu¬ 
lation of the perturbation analysis that could serve as a 
reference for a more complete future study. In particu¬ 
lar, we incorporate the equilibrium inflow speed uor into 
the linearization scheme. We also discuss the expected 
effect of including the equilibrium magnetic tension term 
(oc BzdBr/dz), another key component of the wind¬ 
driving disk model, and show that it could lead to ad¬ 
ditional asymmetric resonances and thereby potentially 
modify the differential torque. However, our applications 
concentrate on the magnetic terms that are relevant to 
vertical angular momentum transport, Bz and dBp,jdz^ 
which we first consider separately and then joi ntly. Our 
semian alytic scheme is patterned on the work of |Terquem| 
(20031, who focused on the case of a pure-H,^ field in 2D. 
(We consider this case and its generalization to 3D in 
Paper I.) For any given field configuration, we identify 
the resonances and the turning points of the governing 
second-order differential equation and locate the wave 
propagation regions where the bulk of the interaction 
between the planet and the disk takes place. In calcu¬ 
lating the torque, we take account of both the point-like 
contribution from the resonance location and the contri¬ 
bution from the surrounding wave-propagation region. 
For a magnetized disk, the relevant waves are SMS and 
Alfven (associated with the magnetic and Alfven reso¬ 
nances, respectively) as well as FMS. The latter gener¬ 
alize the acoustic waves that propagate away from the 
effective Lindblad resonances (which are, in fact, turning 
points) in an HD disk. 

The pure-i ?2 case was previously studied b ylMuto et al.| 
(20081. However, they employed a shearing-sheet formu- 
lation, which does not enable an evaluation of the differ¬ 
ential torque. In addition, they modeled the 3D interac¬ 
tion by using the WKB approximation, which prevented 
them from identifying the full set of turning points for 
this case, and by considering the limit <C 1 , which 
corresp onds to an unrealist ically strong field. Thus, even 
though Muto et al. (2008) identified the relevant waves 
and resonances for this problem, their 3D results are of 
limited applicability and cannot always be directly com¬ 
pared with those obtained using our more general formu¬ 
lation. In the 2D limit we confirmed these authors’ find¬ 
ing that the torque is reduced in comparison with the HD 
regime as the field strength increases, which would have 
the effect of slowing down inward migration. We pointed 
out that, in contrast to the case of a pure-Hta field in 
2D, for which Terquem (2003) and Fromang et al. (2005) 
showed that an outward decrease in the field amplitude 
would tend to reduce the inward drift and might even 
reverse it, for the pure-i?^ case such a decrease would 
have the opposite effect, acting to speed up inward mi- 
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gration. In 3D we found that the contributions of the 
kz ^ 0 modes roughly cancel each other out for typical 
wind-driving disk parameters, so that the total torque 
is comparable to that in the 2D limit, again implying a 
reduced, but still inward-directed, migration. 

We found that a stand-alone dB^/dz field configura¬ 
tion has little effect on planet migration but that the 
torque from a disk in which both Bz and dB^p/dz are 
nonzero is dramatically increased in comparison with the 
pure-U^ case. In the 2D limit, this increase is associated 
in part with the appearance of a wave propagation re¬ 
gion much closer to the planet, where the planet’s grav¬ 
itational influence is correspondingly greater. Thus, for 
a vertical field parameterized by j3z = 0.8, the (dimen¬ 
sionless) total torque increases from ~ 700 in the pure- 
Bz case (as compared with ~ 1300 in an HD disk) to 
^ 3400 when a comparatively weak B^p gradient term 
(corresponding to a surface field \B^ps\ = 0.13B^q) is 
added. In contrast with the pure-B^ case, the net 2D 
torque appears to increase with the field strength when 
the gradient term is added. This suggests that inward 
migration, which speeds up for a uniform field distri¬ 
bution, may be reduced if the field amplitude decreases 
sufficiently fast with radius, although we concluded that 
it is unlikely to be reversed given that in this case the 
contributions to the torque from the inner wave prop¬ 
agation region have the same sign on the two sides of 
the planet. We did not explicitly pursue this possibility 
since we found that the two sides of the planet also con¬ 
tribute with the same sign for the magnetic and Alfven 
resonance regions of the kz ^ 0 modes, where the bulk 
of the torque is produced, which makes the sign of the 
total torque insensitive to any spatial variation of the 
field. For the fcz ^ 0 modes we further found that, for 
high values of the azimuthal mode number m, the MR 
regions become dominant and contribute a strong posi¬ 
tive torque (which induces inward migration). By far the 
largest contribution is produced by the lowest vertical 
mode {kzh = 1.56), for which we obtained a cumulative 
torque of ^ 230,000. We thus infer inward migration at 
a rate that is > 200 times faster than in an HD disk. 
We interpret this result as a manifestation of the ability 
of the planet to plug into the efficient vertical magnetic- 


angular-momentum transport mechanism that operates 
in a disk with this field configuration. This behavior is 
fundamentally different from that of planets in the stan¬ 
dard Type-I migration scenario, which does not involve 
direct coupling to the disk’s underlying angular momen¬ 
tum transport mechanism. 

A more precise determination of the torque in a wind¬ 
driving disk must await a full, self-consistent treatment 
within the framework of 3D, nonideal MHD. In particu¬ 
lar, a numerical simulation that incorporates the disk’s 
magnetic diffusivity could follow the evolution of a planet 
in a quasi-equilibrium wind-driving disk and determine 
not only the implications of vertical angular momentum 
transport but also those of other relevant effects that 
were neglected in the current treatment, such as the mag¬ 
netic tension force and the background accretion veloc¬ 
ity. Furthermore, resistive and viscous dissipation effects 
that become important on small spatial scales could in 
practice limit the effect of the high-m azimuthal modes 
for kzh = 1.56 that, according to the ideal-MHD linear 
analysis, contribute most of the positive torque on the 
disk. It would also be necessary to determine whether the 
strong field-matter coupling condition (Elsasser number 
A ^ I) that underlies our formulation indeed applies at 
the midplane of the planet formation region of real pro¬ 
toplanetary disks. In fact, even if the main conclusion of 
our study is corroborated by a more detailed investiga¬ 
tion, it would be possible for a small planet located in a 
wind-driving region of such a disk to avoid rapid inward 
migration if the above condition is not satisfied at that 
location. 
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Gomputing Gluster at the Research Computing Center 
of the University of Chicago. 


APPENDIX 


A. WAVE PROPAGATION WHEN K IS COMPLEX 


When the coefficient K given by Equation (22) is complex (i.e., if it can be written in the form K = a + ib, where a 
and b are real numbers and 5^0), the solution always has a wave-like component (i.e., it can be written in the form 
C exp [(c -I- id)r], where C, c, and d are real numbers, and d ^ 0). Now, c < 0 corresponds to wave damping over a 
characteristic (e-folding) length ~ I/|c|. A rough criterion for wave propagation in the presence of damping is tha t 
this length be larger than the size of the wave-format ion (resonance) region, which is ^ I/|d| (e.g., Artymowicz|I993 ); 
in other words, a wave propagation region is characterized by I- One can solve for c and d in terms of a and 

b by plugging the adopted form of the solution into Equation (21). This yields 


^ ~ dP- 


c" {-s/ap -I- 52 _ aY 


52 


(AI) 


Equation (Al) can be used to determine the wave propagation regions for given values of a and b that define the 
coefficient K . One can gain insight into this question by considering the behavior of Y in the limits of small and large 
X = \? jaP. When A ^ I, the result depends on the sign of a. For a > 0, U —i.e., F ^ I, implying wave 
propagation. This is consistent with the standard hydrodynamical result that corresponds to 9{B} = 0 (or A = 0), 
where wave propagation is inferred to occur when SfijAT} > 0. On the other hand, when A ^ I but a < 0, F —^ 
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Fig. 9. — Real and imaginary parts of the coefficient K (Eq uatio n ( |22[ |) and the ratio Y = jd? (Equation ||A1[ |) for the 2D limit of the 
i?z 0, dBipJdz 7 ^ 0 field configuration considered in Section |4.4| The regions of wave propagation are identified by the conditions T ~ 1, 
|St{R'}| >> ISRIif}! and Y ~ 0, SRliC} > 0. The evanescence regions correspond to the locations where Y 1. These results are used in 
the construction of the schematic shown in the top panel of Figure]^ 


is 3> 1, indicating a region of evanescence (again in accord with the behavior of a hydrodynamical disk in the limit 

X = 0). 

When X 3> 1, Equation (All implies F —>■ 1, remaining < 1 for a > 0 and > 1 for a < 0. In this case the wave 
propagation criterion formulated above is only marginally satisfied. This case is realized in the vicinity of the planet 
in the 2D limit {kz = 0) of the Bz yf 0, dB^p/dz 0 field configuration considered in Section 4.4 as shown in FigureJ^ 


The profile of the perturbed enthalpy for this model indeed hints at wave-like behavior in that region (see Figure I 
and we have therefore labeled it as a wave-propagation zone in the schematic presented in the top panel Figured^ 
We have delimited the extents of the wave-propagation and evanescence regions in this schematic based on the results 
shown in Figure]^ 


B. DOMINANT TORQUE TERMS FOR THE “COMBINED” FIELD CONFIGURATION 

We wish to determine whether a planet located in a disk with a combined Bz + dB^^/dz field configuration indeed 
taps into the vertical magnetic transport of angular momentum (represented by the oc BzOB^/dz term in the angular 
momentum conservation equation). To this end, we examine the dominant terms in the expr essi on for the perturbed 
magnetic torque density acting on the disk, (r — rp)F^. From the ip component of Equation (|14[) we have 


piF'^ = B'^dBo^/dz -f i{kzB'p, - kpB'^)Boz 


(Bl) 


where k^ = m/r. In view of our finding in Section 


4.4 


that the behavior of the disk is different for the kz = 0 and 
kz ^ 0 vertical modes, we consider these two cases separately. 

For kz yf 0, the bottom entry of Table [T] highlights the fact that high azimuthal mode numbers m dominate the 
torque. It is also clear from this entry that the regio n around the magnetic resonance provides the main co ntrib ution 


to the torque. Therefore, we consider Equation (Bl) in the MR regio n an d in the limit of a large k^p. Figure 10 shows 


the behavior of the three terms on the right-hand side of Equation (Bl| in the inner MR region for the parameters 


that correspond to the aforementioned table entry: it is seen that the last term, which contains an explicit dependence 
1 kp, clearly dominates. 

The z component of the perturbed induction equation (Equation §) gives 


ik^raB'^ = 


Substituting the perturbed continuity equation, 


dr 


Boz - ikpV^ 




Oz ■ 


r I -r dv’ 

ikpva— -I--h 

Po \ r dr 


+ ik^v'^p F ikzv'^ = 0 , 


(B2) 


(B3) 
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Fig. 10.— Behavior of the (real) contributions to the azimuthal component of the perturbed Lorentz force density (Equation ( |B1[ |) in the 
vicinity of the inner magnetic resonance for the model parameters that correspond to the bottom entry of Table [b The term oc Ic^B'^Bqz 
is seen to dominate at and around the resonance. 


into Equation (B2), we arrive at the following form for the perturbed vertical field component: 

kz 


Po 


k^pva 


v'Bt 


Oz ■ 


(B4) 


Since we are focusing on the regions around the magnetic resonances (which lie close to the planet), a is generally 
small, so even for large values of one can expect the second term on the right-hand side of Equation (B4) to be the 
principal component of B'^. We have verified numerically that this term is indeed the largest factor in the expression for 

B'^ (although the contribution of the hrst term remains nonnegligible). In our analytic formulation, p /po (x W is given 
by a complicated expression in which seve ral t erms depend on the product B^^dBi^^pldz. However, the dependence of 
the dominant (second) term in Equation (B4) on this product is more clear-cut. The perturbed vertical velocity, u(,, 
is induced primarily by the vertical component of the Lorentz force, which for this held conhguration is proportional 
to {B'^dB^^/dz)/ pq. Thus, the leading term that affects the angular momentum transfer from the planet to the disk 
in this case is (substituting the above dependencies into Equation (Bl)): 




kzB'^B^z 

Po 




dB, 


Oz- 


Tty 


(B5) 


It is evident from Equation (B5) that the vertical transport of angular momentum by the held plays a major role in 
shaping the 3D behavior of a disk with this hel d con hguration. 

When kz = 0, the middle term of Equation (Bl) vanishes, and we are again left with terms that depend on the 
perturbed vertical held (which in this case is given by B'^ = Boz:{p'/po)). Since the interesting behavior of the 2D 
case illust rated in Figure]^ is found to be at low m (where the inner wave propagation region is largest), we look at 
Equation (Bl) at small In this limit the B'^dB^^/dz term is dominant. Accordingly, oc {p'jpt^B^^dB^^Idz, 


demonstrating that the vertical magnetic transport of angular momentum plays a role in the 2D limit as well. 

Our hnding that the dominant contributions to the torque for both the kz = 0 and kz 0 modes have the same sign 
on the two sides of the planet is consistent with the picture that in the “combined” field case the planet loses angular 
momentum primarily through the large-scale magnetic field. The planet couples to the held indirectly, through the 
density perturbations that are induced in the gas by its gravitational potential. The perturbed held exerts a back 
torque that is transmitted to the planet (again, through its gravitational interaction with the gas) even as the held 
transports the planet’s angular momentum to the disk surfaces, where it can be deposited in a centrifugally driven 
wind (or, alternatively, in torsional Alfven waves that propagate into the ambient interstellar medium). In our formal 
treatment, the angular momentum transport is described through the launching of MHD waves that propagate both 
radially and vertically into the surrounding gas. However, the dominant physical transport is ehected by the large-scale 
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field that carries the angular momentum in the vertical direction. And while the analysis of the wave propagation 
regions in the disk is useful for identifying the locations where the planet interacts most strongly with its surroundings, 
the standard picture in which the torque exerted on the planet depends on the planet’s angular velocity relative to 
the gas (and therefore changes sign between the interior and exterior disk regions) does not apply in this case. In 
the standard picture, Type I migration does not depend on the underlying (“viscous”) angular momentum transport 
mechanism in the disk. Th e effective disk kinematic viscosity v is commonly represented in the form v = 

(Shakura & Sunyaev 1973), with the parameter a typically taken to be ^ 1, and its relative unimportance in the 
Type-1 migration process is due to the fact that the viscous transport term is a factor ^ a smaller than the other terms 
in the perturbed angular momentum conservation equation. By contrast, for disks in which a large-scale magnetic 
field dominates the angular momentum transport, the effective value of a is 1, so in this case the disk’s angular 
momentum transport mechanism can directly affect the planet’s migration. 
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